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Abstract. Numerical algorithms for solving problems of mathematical physics 
on modern parallel computers employ various domain decomposition tech- 
niques. Domain decomposition schemes are developed here to solve numer- 
ically initial/boundary value problems for the Stokes system of equations in 
the primitive variables pressure-velocity. Unconditionally stable schemes of 
domain decomposition are based on the partition of unit for a computational 
domain and the corresponding Hilbert spaces of grid functions. 



1. Introduction 

In computational fluid dynamics [Jl llOj there are employed numerical algo- 
rithms based on using the primitive variables pressure-velocity. The main difficul- 
ties in this approach are connected with the calculation of the pressure. In studying 
transient problems the corresponding elliptic Neumann problem for the pressure is 
derived as the result of employment of one or another scheme of splitting with 
respect to physical processes [3l [4] . 

Domain decomposition methods are used for the numerical solution of bound- 
ary value problems for partial differential equations on parallel computers. They 
are in most common use for stationary problems O 111] . Computational algo- 
rithms with and without overlapping of subdomains are employed in this case in 
synchronous (sequential) and asynchronous (parallel) algorithms. 

For transient problems it seems to be more suitable to utilize iteration-free 
variants of domain decomposition techniques [6l [8] which are best suited to pe- 
culiarities of a problem (evolution in time). In these regionally- additive schemes 
a transition to a new time level is performed via solving problems in particular 
subdomains. 
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The regionally-additive schemes for the Navier-Stokes equations in the primi- 
tive variables are discussed in [2]. In simulation of incompressible flows an ellip- 
tic problem for the pressure can be changed to separate elliptic problems for the 
pressure in particular subdomains. Therefore, it is possible to construct iteration- 
free regionally-additive schemes for the Navier-Stokes equations. In this paper we 
propose a general approach to construct domain decomposition schemes for time- 
dependent systems of equations. Using the partition of unit for a computational 
domain and the corresponding Hilbert spaces of grid functions we perform a tran- 
sition to finding the individual components of the solution in the subdomains. The 
unsteady Stokes equations for an incompressible fiuid is considered as a typical 
problem. 



2. Stokes equations 

Assume that the linear approximation is valid to describe a flow of an incom- 
pressible viscous fluid. In a region with solid boundaries we can write equations 
of motion and continuity in the primitive variables pressure, velocity as follows 

(2.1) ^+gradp-j. Au = f(x,i), 

(2.2) divu^O, xefJ, 0<i<T. 

Here u is the velocity, p is the pressure, v is the kinematic viscosity and A = 
div grad is the Laplace operator. Equations (|2.ip . (|2.2p are supplemented with the 
following condition for the single- valued evaluation of the pressure 

(2.3) j p(x, t)dx = 0, < t s; T. 

n 

No-slip, no-permeability conditions are specified on solid boundaries 

(2.4) u(x, t)=0, X e < i sC T. 
Some initial condition is also given 

(2.5) u(x,0)=v(x), xerj. 



Let us rewrite problem (j2.ip - (j2.5p in an operator formulation. On the set of 
functions satisfying (|2.3p . [^^ . we have the Cauchy problem 

(2.6) ^ + Au + Bp = {, 

at 

(2.7) S*u = 0, 0<t!^T, 

(2.8) u(0) = V. 
For these operators in the space L2(r2) we have 

A^A*^S£, S = 5{n) > 0, 
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where £ is the unit (identity) operator. Adjointness of operators B — grad {B : 
L2{n) L2in)) and B* = - div (B* : L2{n) L2{n)) fohows from 

/ ugradp(ix+ / divupdx^O. 

For problem (j2.6p - (l2.8p the fohowing simple a priori estimate is valid 

(2.9) ||u(Of ^||vf + l^Vwfrf^, 

where || • || is the norm in L2(il). Estimate p.9p will be for us a reference point 
when considering discrete problems. 

3. Approximation in space 

In this work the main attention is paid to computational algorithms for the 
transition to a new time level, i.e. approximation in time. To construct discretiza- 
tion in time, operator-splitting schemes are used that allows to formulate a problem 
for the pressure in the most natural way. The problem of approximation in space 
is solved in the standard manner. 

There are employed various types of grids: the non-staggered (collocated) grid 
where both the pressure and velocity components are referred to the same points; 
next, partially staggered (ALE- type) grid where the pressure is referred to an in- 
dividual grid shifted in all space directions on a one-half step from the basic grid 
where all velocity components are defined; and finally, the staggered (MAC-type) 
grid where the pressure is defined at the center points of grid cells whereas the 
velocity components are referred to the corresponding faces of the cell. 

For simplicity we consider here uniform rectangular non-staggered grids. Prob- 
lem (|2.ip - (|2.5p is solved in a rectangle 

= { X I X = {xi,X2), < Xa < la, a = 1, 2}. 

The approximate solution is calculated at the points of a uniform rectangular grid 
in J7: 

Ol = {x I X = (a::i,X2), Xa ^ iaha, ia ^0,l,...,Na, Naha^la} 

and let uj be the set of internal nodes {lj = uj U duj). 

For vector grid functions u(x) = 0, x e dijj we define a Hilbert space H = L2(a;) 
with the scalar product and norm 

(U, V) = 5] U(x)v(x)/ii/l2, ||U|| = (U, u)l/2. 

The grid operator A is taken in the form A = —vl^h, where A/i is the grid Laplace 
operator: 

^hU = -i(2/(xi -I- hi,X2) - 2y{xi,X2) + y{xi - hi,X2)) 
"-1 

- -72 (y (2^1, 3^2 + h2) - 2y{xi,X2) + y{xi,X2 - h2)) ■ 
"2 
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In H the operator A is selfadjoint and positive definite: 



2 . , 

(3.1) A = -i.A, =A* ^i^ShE, 5, = sin^^. 

a=l " " 

The pressure gradient is approximated by directed differences with an error 0(h). 
We set B ~ grad^ at 

(3.2) Bp = {{Bp)i,iBp)2}, xGc^, 
where 

{Bp)i = —{pixi + hi,X2) -p{xi,X2)), 

hi 

{Bp)2 = ■j-{p{xi,X2 + ft.2) -p(xi,X2)), 

ll'2 

The set of points for the pressure evaluation is denoted as Wp (wp C ui). For the 
grid divergence operator B* = — div/i we have 

(3.3) B*VL = --^{ui{xi,X2) - Ui{xi - hi,X2)) + 

III 

~J-{U2{X1,X2) - U2(xi,X2 - ^-2)), X G LOp. 
h2 

The adjointness property of the grid gradient and divergence operators is a conse- 
quence of the discrete equation 

Sp(x)u(x)(x)/ii/i2 + p(x)B*u(x)(x)/ii/i2 = 0, 

which takes place on the set of vector grid functions u(x) = 0, x G duj. 

In view of p.l[) - (l3.3p after the spatial approximation of problem (|2.6p - (|2.8p we 
obtain the following problem 

(3.4) ^+Au + Bp = {, 

(3.5) B*u = 0, < < ^ T, 

(3.6) u(0) = V. 

For the solution of problem p.4p - (|3.6p a priori estimate (12. 9p holds, where now || ■ || 
is the norm in H = L2(a;). 

4. Domain decomposition 

Let be a combination of p particular subdomains 

17 = f^i U 1^2 U ... U fl,n ■ 

Particular subdomains can overlap one onto another. We shall construct the 
schemes of decomposition where the solution at the new time level for the initial 
problem is reduced to the sequential solution of problems in particular subdomains. 
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Let us define functions for domain D, 

(4.1) = { ^o!'' x^^^a!r a^l,2,...,m. 

Generally, see for example, [6], [8], domain decomposition schemes for unsteady 
problems are based on the partition of unit for the region fl, where 

Q = l 

It is more convenient to use a somewhat different partition where 

rn 

(4.2) ^772(^)^1^ 

For the decomposition of computational domain (14. ip . (|4.2p we consider the 
following additive representation of the identity operator E in H = L2(w): 



(4.3) E^^xl, Xa^Va{^)E, X6W, a = l,2,...,m. 

a=l 

Taking into account (14. 3p we have 

m 

(4.4) u = ^Uq, Ua = XaU, a = l,2,...,m. 

a=l 

To formulate an appropriate system of equations for determining components 
of the solution Uq , a = 1,2,..., 171^ we multiply both, sides of equation (|3.4p by x.a- 
This gives 

(4.5) -^+XaA^Xf!Ul3+BaP^ia: 

13=1 

where 

Ba^XaB, fa = Xaf, a = l,2,...,m. 
Taking into account that 

equation fj3.5p in the new notation is written as 

(4.6) B*u„ =0, < t T. 

The system of equations (|4.5p . (14. 6p is supplemented by the initial conditions 

(4.7) Uq(0)=Vc, Va = XaV, a =1,2,..., TO. 

For C/ = {ui, U2, u,„}, we define in the space H™ the norm and inner product 
as follows 



(c/,y)^ = ^(u„,v„), ||c/||^ = (c/,c/) 



1/2 
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Let US multiply scalar ly individual equations (|4.5p by Uq,, a = 1,2,..., to and add 
them together. Next, multiply equation (j4.6p by p. Taking into account (|4.4L we 
obtain 

a— 1 Q— 1 

This implies the a priori estimate 

(4.8) ||;7||^scexp(i)||y||^+ f eMt-e)\\F{e)\\lde 

for problem with = {fi, fa, £,„}. 

5. Splitting scheme 

In the construction of domain decomposition schemes we shall proceed from 
the scheme of splitting with respect to physical processes for the Cauchy (|2.6p - 
(|2.8p . We shall use a simple additive scheme componentwise splitting [H [9]. Let 
u" be the difference solution at the time moment = nr, where t ~ T/N > is 
the time-step. Let us separate out a particular stage connected with the pressure 
impact [3l |4j . Thus, in the first stage we have 

,,n+l/2 _ „ji 

(5.1) — + Au"+i/2 = f"+i/2. 

T 

The pressure gradient is treated only in the second stage: 

„n+l _ „n+l/2 

(5.2) + = 0, 

T 

(5.3) B*u"+i = 0. 

Implementation of (|5.2p . (j5.3p consists of two steps. In the first step we solve the 
following problem for the pressure 

T 

whereas in the second one we update the velocity: 
Multiplying (|5.ip by u"+^/^, we obtain 

||U"+1/2||2<; ||u«||2 + ^||p+l/2||2^ 
Oh 

Similarly, from (j5.2p taking into account (j5.3l) we have 

||^„+l||2 ^ ||un+l/2||2^ 

Thus, we obtain the grid analog of estimate (|2.9p 

||u"+l||2^||u"f + -^||r+l/2||2_ 

for difference scheme (|5.ip - (|5.3p . 
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For simplicity, we shall construct splitting schemes for problem (|4.5p - (|4.7p by 
analogy with splitting scheme (|5.ip - (|5.3p . The first half-step (viscous dissipation) 
is associated with the solution of equations 



Xa^^X/3U^ =fa, a = l,2,...,m, r<<^f"+i/^ 



3=1 



Taking into account the fact that the numerical solution is implemented via solving 
individual problems in the subdomains, the transition from time level i" to level 
^g^^ realized as follows: 

ji+l/4 _ „ g-i , 

(5.4) + x.A Y: XP^T" + ixo^X.<+V4 ^ 

a — 1,2, m, 

n + l/2 Ti+1/4 , m 

(5.5) ]^ + ^XaAxa.K+'/'+XaA ^ X/^u^+^Z' = 0, 

/3=a+l 

a = 1, 2, m. 

In view of (j5.4p . (j5.5l) in each subdomain JIq, a = l,2,...,m we must invert grid 
selfadjoint elliptic operator 

Da = E+ ^XaAXa 

for finding u"^'^^^ {a — l,2,...,m) and Uq^^^^ {a — 171,171 — 1,...,1). In this case, 
outside subdomains fla, = 1, 2, m there are used explicit calculations. 

Stability of scheme (|5.4p . (|5.5p will be investigated in H™. Consider the oper- 
ator 

(5.6) A^lAafi}, Aap = XaAxi3, a, /3 = 1, 2, 771. 

Taking into account ([XT]) , (g^), we have A = A* > in H". Scheme EH), (E3 is 
based on using the triangular splitting 

(5.7) A = Ai+A2, Ai=A;. 
Using notation (|5.6p . (|5.7p we rewrite [5^ . (|5.5I) in the form 

Tjn+l/i _ jjn 

(5.8) — + Ai;7"+i/* = 

T 

(5.9) + A2C/"+i/2 = 0. 

T 

Taking into account that Aq ^0, a = 1, 2 in H™, for (15. 9p we immediately have 

(5.10) 

Muhiplying ([5?S)) by [/"+i/4^ .^^g obtain 
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For the last term on the right hand side we use the estimate 

2t{F-+^'\ t/"+V4)^ (1 _ exp(-r)) ||L/"+i/4||2 ^ WF-^''^. 

1 - exp(-T)) 

This leads to the estimate 

(5.11) \\U-^"'\\l < exp(r) ||C/"|1^ + r||F"+V2 1|2^. 

The second half-step results from the pressure and is connected with the system 
of equations 

—j—+BaP = 0, a = l,2,...,m, 
at 

m 

Q = l 

Approximation in time for such systems were considered in |12j . We shall use the 
additive scheme 

(5.12) u«+l/2+/3/2™ ^ ^n+l/2+(/3-l)/2™^ ^ /? = 1 , 2, . . . , m, 

n-|-l/2+Q/2m n+l/2+(Q- l)/2m 

(5.13) ^^^^ ^-^^^^ + S^p"+l/2+a/2m ^ g, 

T 

(5.14) B*<+i/2+"/2m ^ 0, a = l,2,...,m. 

The implementation of additive scheme (|5.12p - (|5.14p is conducted by analogy with 
scheme ((O)) . 

For ([ET^ - ([ETi|) we have 

IK+^II ^ IK+i/^ll, a = l,2,...,m 

and therefore 
(5.15) 

Taking into account (j5.10p , (j5.1ip and (I5.15P we obtain the desired stability estimate 
of the additive operator-difference scheme (15.41) . (15. 5p . (|5.12l) - (|5.14p 

(5.16) < exp(r) ||C/"||?„ + 

Estimate (|5.16p is the grid analog of (|4.8I) for the differential problem. 
This allows us to formulate the following main result. 

Theorem 5.1. The additive scheme of domain decomposition |5.4[ ), i5.5]) . liS.l^) - 
ji5.14\ l is unconditionally stable and estimate i5.16]) holds for the numerical solution. 
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